home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / FFTSING.ARJ / VERSING.C < prev    next >
C/C++ Source or Header  |  1992-04-15  |  5KB  |  197 lines
  1. /*
  2.  *    versing.c   Verifies sing.c  
  3.  *
  4.  *    Generates a square pulse of width width and total length of
  5.  *   2*half_order. Computes the discrete Fourier transform and compares
  6.  *   it to the theoretical transform.
  7.  *   Obtains the inverse transform and compares it to the original pulse.
  8.  */
  9.  
  10.                 /* Version 2.0  April 1992 */
  11.  
  12. /**************************************************************************
  13.  
  14.     Javier Soley, Ph. D,   FJSOLEY @UCRVM2.BITNET
  15.     Escuela de Física y Centro de Investigaciones Geofísicas
  16.     Universidad de Costa Rica
  17.  
  18. ***************************************************************************/
  19.  
  20.  
  21. /* Includes */
  22.  
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <alloc.h>
  27. #include <string.h>
  28. #include <dos.h>
  29. #include <time.h>
  30.  
  31. /* Defines */
  32.  
  33. #define    TWO_PI    ((double)2.0 * M_PI)
  34.  
  35. /* Globals */
  36.  
  37. unsigned long  half_order,  mem_avail, max_size, width;
  38. double huge *hreal, huge *himag;                 
  39.  /* use to reference long series without
  40.                           wrap around */
  41.  
  42. /* Prototypes and forward declarations */
  43.  
  44.  
  45. double amp_rect_pulse (unsigned long, unsigned long, unsigned long);
  46. double pha_rect_pulse (unsigned long, unsigned long, unsigned long);
  47.  
  48.  
  49.  
  50. /* The program */
  51.  
  52. void  main(void)
  53. {
  54.     double far  *real;
  55.     double far  *imag;   /* returned by memory allocation functions */
  56.     unsigned long counter, j,  r_index, i_index;
  57.     long start;
  58.     double tmag, tpha;  /* magnitud and phase of theoretical FFT */
  59.     time_t start_time, finish;
  60.  
  61.     mem_avail = farcoreleft();
  62.     printf("%lu bytes free in far heap\n", mem_avail);
  63.     max_size =  (mem_avail-2) / sizeof(double) ;
  64.     printf( "Maximum length of data is %lu \n", max_size);
  65.     printf( "Maximum half-length is %lu  \n", max_size /2);
  66.  
  67.     printf("Half-length of rectangular pulse ?  ");
  68.     cscanf("%lu", &half_order);
  69.     if ( 2*half_order  > max_size) {
  70.         printf("\nLength of series exceeds far heap\n");
  71.         exit(2);
  72.     }
  73.  
  74.     printf("\nWidth of rectangular pulse ?  "); cscanf("%lu", &width);
  75.     if ( width  >= 2*half_order) {
  76.         printf("\nWidth of pulse exceeds length of pulse\n");
  77.         exit(3);
  78.     }
  79.  
  80.     
  81.     if ((real =  farcalloc(half_order +1, sizeof(double))) == NULL) {
  82.         printf("\nAllocation error with real part\n");
  83.         exit(4);
  84.     }
  85.  
  86.     if ((imag =  farcalloc(half_order + 1, sizeof(double))) == NULL) {
  87.         printf("\nAllocation error with imaginary part\n");
  88.         exit(5);
  89.     }
  90.     hreal = real; himag = imag;  
  91.     printf("\nThe real e imaginary parts start at %Fp and %Fp\n",
  92.         hreal, himag);
  93.  
  94. /* Generate the rectangular pulse */
  95.  
  96.     r_index = 0; i_index =0;  counter =0;
  97.     while ( counter < width ) {
  98.         hreal[r_index] = (double) 1.0;
  99.         if (counter == width -1) break;
  100.         r_index++; counter ++;
  101.         himag[i_index] = (double) 1.0;
  102.         i_index++; counter++;
  103.     }
  104.     r_index++;
  105.     while( r_index <= half_order ) {
  106.         hreal[r_index] = (double) 0; r_index++;
  107.     }
  108.     i_index++;
  109.     while( i_index <= half_order ) {
  110.         hreal[i_index] = (double) 0; i_index++;
  111.     }
  112.  
  113.         /* Now call Singleton's FFT procedures */
  114.     start_time = time(NULL);
  115.     sing(hreal, himag, half_order, half_order, half_order, -1);
  116.     realtr(hreal, himag, half_order, -1);
  117.     finish = time(NULL);
  118.     printf("\nIt took aproximately %f  seconds\n",
  119.      difftime(finish, start_time)); 
  120.                 /* Whistle when done */
  121.     sound(440); delay(4000); nosound();
  122.  
  123.     do {
  124.         printf("\n\nNegative number to exit");
  125.         printf("\nVerify results starting where ? [0..%lu ] ",
  126.               half_order-19);
  127.         cscanf("%lu", &start); 
  128.         if( start <0 ) break;
  129.         clrscr();
  130.         printf("          Real theory     Real calculated    Imag. theory     Imag.calculated\n");
  131.  
  132.             for ( j=start; j < start + 20; j++) {
  133.                 tmag = amp_rect_pulse(j, width, 2*half_order);
  134.                 tpha = pha_rect_pulse(j, width, 2*half_order);
  135.                 printf( "\n%5lu  %16.9e  %16.9e  %16.9e  %16.9e", j,
  136.                 tmag*cos(tpha), hreal[j]/(4*half_order),
  137.                 tmag*sin(tpha), himag[j]/(4*half_order));
  138.             }
  139.        } while (1);
  140.  
  141.         /* Now call the inverse transform */
  142.     printf("\nCalculating the inverse transform\n");
  143.     start_time = time(NULL);
  144.     realtr(hreal, himag, half_order, 1);
  145.     sing(hreal, himag, half_order, half_order, half_order, 1);
  146.     
  147.     finish = time(NULL);
  148.     printf("\nIt took aproximately %f  seconds\n",
  149.      difftime(finish, start_time)); 
  150.                 /* Whistle when done */
  151.     sound(440); delay(4000); nosound();
  152.  
  153.     do {
  154.         printf("\n\nNegative number to exit");
  155.         printf("\nVerify results starting where ? [0..%lu ] ",
  156.               half_order-20);
  157.         cscanf("%lu", &start); 
  158.         if( start <0 ) break;
  159.         clrscr();
  160.         printf("            Even index       Odd index \n");
  161.  
  162.             for ( j=start; j < start + 20; j++) {
  163.                 printf( "\n%5lu  %16.9e  %16.9e ", j,
  164.                 hreal[j]/(4*half_order),himag[j]/(4*half_order));
  165.             }
  166.        } while (1);       
  167.                     
  168. }
  169.  
  170.  
  171.  
  172.  
  173. /* Calculates the amplitude spectra of a rectangular pulse */
  174.  
  175. double amp_rect_pulse (i, w, n)
  176. unsigned long i, w, n;  /* frequency index, width of pulse, order of fft)*/
  177. {    double angle, amp;
  178.     angle =     M_PI*i/n;
  179.     if ( i==0 || i==n) amp = (double) w / (double)n;
  180.     else     {
  181.     amp = sin(w*angle)/sin(angle);
  182.     amp /= n;
  183.     }
  184.     return amp;
  185. }
  186. /* Calculates the phase spectra of a rectangular pulse */
  187.  
  188. double pha_rect_pulse (i, w, n)
  189. unsigned long i, w, n;  /* frequency index, width of pulse, order of fft)*/
  190. {    double phase;
  191.     phase = (double)M_PI*(double)(w-1)* (double)i/(double)n;
  192.     return -phase;
  193. }
  194.     
  195.  
  196. /* EOF */
  197.